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I. INTRODUCTION 

Several outstanding questions in cosmology, such as 
the current accelerated expansion of the Universe [1] or 
the presence and origin of dark matter [2], have motived 
many efforts to expand upon the conventional formalism 
of general relativity (GR), as the latter cannot provide a 
direct explanation of such phenomena. A great number 
of attempts have been made to phenomenologically study 
putative modifications of GR, ranging from the inclusion 
of nontrivial curvature terms in the Einstein-Hilbert ac¬ 
tion [3, 4] to the addition of suitable scalar [5] and vector 
fields [6-9], among others [10]. 

Amongst the latter, the bumblebee model was initially 
put forwarded in 1989 [11] as a toy model for spontaneous 
Lorentz symmetry breaking using a vector field; being 
an extension of GR, it takes the usual Einstein-Hilbert 
action and expands it with a vector field and a potential. 
It provided a simple and more tractable scenario than the 
more general Standard Model extension, a framework for 
studying Lorentz-breaking terms added to the Standard 
model of interactions and general relativity [12]; several 
aspects of the impact of the SME on the gravitational 
sector have been addressed e.g. in Refs. [13-18]. 

The bumblebee model can be considered as a particu¬ 
lar case of vector-tensor models [6, 7], which also include 
the more well-known Einstein-aether models [8]; in the 
latter, the potential of the vector field resorts to a La¬ 
grange multiplier that constrains its norm, and exhibits 
couplings between the vector field and the Ricci tensor 
and scalar. Distinctively, the bumblebee model assumes 
only that the potential has a nonvanishing vacuum ex¬ 
pectation value (YEV), but restricts the coupling to just 
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the Ricci tensor, as depicted in the action functional (1). 

The bumblebee model incorporates a mechanism of 
spontaneous breaking of the Lorentz symmetry, which 
is compatible with general Riemann-Cartan geometries 
[9, 12, 19-22]. This mechanism is inspired by the Higgs 
mechanism of the Standard model of fundamental par¬ 
ticles and interactions (itself an allusion to the Landau 
theory of phase transitions that first originated in Con¬ 
densed Matter Physics), and serves to both safeguard 
Lorentz symmetry as a symmetry of the action and reg¬ 
ulate the way in which that symmetry is broken. 

The coupling of curvature to the vector field and the 
presence of a potential leads to a higher complexity of 
the equation of motion that cannot be extrapolated from 
approaches that eschew the use of a potential [23, 24] or 
uncouple the vector field from the curvature [25]. 

The aim of this work is to further assess the impact of 
the bumblebee model and explore the effects of a spon¬ 
taneous LSB in a cosmological context — namely its po¬ 
tential as a candidate for dark energy. 

This work is organized as follows: first, the model 
and its governing equations are described; the impact 
of the latter on cosmology is then presented, prompting 
the study of the ensuing dynamical system; finally, the 
obtained results are discussed and conclusions are drawn. 


II. THE BUMBLEBEE MODEL FOR 
SPONTANEOUS LORENTZ SYMMETRY 
BREAKING 

As mentioned in the previous section, the bumblebee 
model extends the standard formalism of general rela¬ 
tivity by allowing a LSB; this is dynamically driven by 
a suitable potential exhibiting a nonvanishing VEV, so 
that the bumblebee vector field B M acquires a specific 
four-dimensional orientation. 

The action functional should capture the most rele- 
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vant features of the bumblebee model, namely the cou¬ 
pling between the bumblebee field and geometry and the 
presence of a general potential. As discussed in detail 
in Ref. [7], the simultaneous inclusion of kinetic terms 
for the vector field of the forms (with F ^ the 

field-strength) and (V^H ^) 2 leads to the appearance of 
ghost degrees of freedom, which vanish if only one of the 
latter is considered. 

Furthermore, an extensive discussion of the symme¬ 
tries of a model with couplings of the form B [1 B V R IIV 
and B^B^R was presented in Ref. [14]: in particular, if 
the potential is quadratic and includes a Lagrange mul¬ 
tiplier A, V(X,B fl ) = A (B^B^ ± 6 2 ), so that variation 
with respect to the former fixes B^B^ = ± 6 2 , then the 
coupling term B^B^R can be gauged away by rescaling 
G. 

This said, the fixed points obtained in this study indi¬ 
cate that the bumblebee field does not rest at its nonvan¬ 
ishing VEV, for a general potential, so that the inclusion 
of both types of couplings should have physical conse¬ 
quences. However, for simplicity we drop the coupling to 
the Ricci scalar and thus consider the model posited in 
Ref. [9], thus obtaining a more computationally tractable 
problem which still retains the main features mentioned 
above; a discussion of the effect of an additional cou¬ 
pling with the Ricci scalar can be found in the Appendix, 
where it is shown that — although the strength of the 
additional coupling could lead to new fixed points — the 
overall dynamical structure of the field equations is not 
substantially modified. 

Given the above, we consider the action functional 


s= / v=9 


2 K 


(R + ^B V R^) - 


4 B ^B, V 


V (B"B m ± b 2 ) + C M 


( 1 ) 


d A x , 


Variation of Eq. (1) with respect to the bumblebee 
field yields its equation of motion, 

= 2 (v'B’' - . (3) 

If the lhs of the equation vanishes, the above results in 
a simple algebraic relation between the bumblebee, its 
potential and the geometry of spacetime. 


III. COSMOLOGY 

In most cosmological studies, the Friedmann- 
Robertson-Walker metric (FRW) metric is assumed, re¬ 
flecting the cosmological principle which posits an ho¬ 
mogeneous and isotropic Universe. However, if Lorentz 
symmetry is spontaneously broken, it is possible that the 
bumblebee field acquires a nonvanishing spatial orienta¬ 
tion, which would dynamically break the aforementioned 
isotropy and require a more evolved geometry. This pos¬ 
sibility (elaborated in the final section) shall not be pur¬ 
sued in this study: instead, one considers the ansatz for 
the bumblebee field, 

B„=(B(i),0) , (4) 

thus upholding the validity of the cosmological principle, 
i.e. maintaining the assumption of large-scale homogene¬ 
ity and isotropy of the Universe. One thus adopts the flat 
FRW metric, as given by the line element, 

ds 2 = — dt 2 + a(t) 2 [dr 2 + r 2 d0 2 + r 2 sin 2 (#)d0 2 ] , (5) 

where a(t ) is the scale factor. 

From Eq. 4, one sees that the field-strength tensor 
vanishes, B^ v = 0 , and the only nontrivial component of 
the bumblebee Eq. (3) is 


where n = 87 rG, £ is a coupling constant (with dimensions 
[£] = M -2 ), Bn is the bumblebee field (with [B^] = M), 
B^ u = d^B u — d v B^ is the field-strength tensor, b 2 = 
b^b 11 = (B^B ^)0 7 ^ 0 is the expectation value for the 
contracted bumblebee vector, V is a potential exhibiting 
a minimum at B^B U ±b 2 = 0, and Cm is the Lagrangian 
density for the matter fields. 

By varying Eq. (1) with respect to the metric, the 
modified Einstein equations are obtained, 
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Hv •> 


( 2 ) 


where V' denotes the derivative of the potential V with 
respect to its argument and is the energy-momentum 
tensor of matter. 



B = 0 


(6) 


which, for a nonvanishing bumblebee field, establishes a 
relation between the dynamics of the potential and the 
scale factor. This is probably the reason that motivated 
the authors of Ref. [7] to describe bumblebee models 
with only a Maxwell kinetic term as “non-dynamical”; 
we maintain, however, that the nonminimal coupling of 
curvature in our case is enough to drive a meaningful 
process, as shall be studied below. 

Considering the matter energy-momentum tensor for 
pressureless dust, 


'CfXU pUfj^Ujy , (7) 

where p is the energy density of matter, u M is the four- 
velocity (with un = ( 1 , 0 ) due to the normalization con¬ 
dition = — 1 ), and using the above expression to¬ 

gether with Eq. (4), the t — t component of the modified 
field Eqs. (2) becomes 

H 2 (l-£B 2 ) = ^(p + V)+ZHBB , ( 8 ) 
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while the diagonal i — i components read 

(ff 2 + 2^ ( 1 -£B 2 ) = (9) 

kV + £ [AHBB + B 2 + BB^j , 

where H = a/a is the Hubble parameter. 

Using the Bianchi identities, one also obtains the fol¬ 
lowing modified equation of conservation of energy, 

p= -3Hp-3^-B{HB + B) + 3^-B 2 , (10) 

Hi a Hi a 

showing that there is an energy exchange between matter 
and the bumblebee field. 

Due to the complexity of the equations derived in this 
section, a simple, closed-form solution cannot be ob¬ 
tained by purely analytical means. Given this difficulty, 
it becomes more enlightening to study instead the full 
dynamical picture presented by the Friedmann and Ray- 
chaudhuri equations under the constraints imposed by 
the bumblebee field equation of motion. 


A. Absence of coupling between the bumblebee 
field and the Ricci tensor 


At this point, it is useful no notice the relevance of the 
coupling strength £: if this is “switched off” by setting 
£ = 0, the modified Einstein Eqs. (2) become 


— T/jlv + 


( 11 ) 


2 V'BpBv - B m BZ - 


V 


\B a pB 


9 hlv 


while the bumblebee equation of motion (3) reads 

= 2 V'B V . (12) 

Although simplified, the above equations in principle still 
allows for dynamical behavior. However, the assumption 
Eq. (4) leads to a vanishing field force B^ v = 0, so that 
the bumblebee Eq. ( 6 ) collapses to 

V'B = 0 , (13) 

forcing the bumblebee field to either vanish or rest at 
one of the extrema of its potential, and keeping it from 
evolving with time. 

Thus, setting £ = 0 and B = B = 0 in Eqs. (8-10) 
yields the simplified set, 

H 2 = 1«(P + V 0 ) , (14) 

H 2 + 2 ^ = kV 0 , (15) 

p = -3 Hp . (16) 

Thus, the bumblebee field contributes solely through the 
value of its constant potential V(B 2 ) = Vq: checking 


the action functional ( 1 ), one sees that it collapses to 
the Einstein-Hilbert form with a Cosmological Constant 
A = kVo- Indeed, Eq. (15) can be directly integrated, 
yielding a de Sitter solution 

a(t) = a 0 e Ho(t - to) . (17) 

with JJq = k.Vq/3 = A/3; Eq. (14) then requires that 
p = 0 , as expected. 

As the next section will show, the existence of a phase 
of exponential acceleration of the Universe is not re¬ 
strained to the uncoupled scenario £ = 0 , but for a non¬ 
vanishing coupling as well. 


B. Generality of a de Sitter solution 


Following the preceding discussion, we consider a so¬ 
lution of the form Eq. (17) (with unknown Ho), so that 
the bumblebee field equation ( 6 ) becomes 

V'B 2 = --H 0 2 B 2 . (18) 

For a nonvanishing bumblebee field, any potential V 
which varies nonlinearly yields a quantity V' that de¬ 
pend on B , and the above implicitly establishes that the 
bumblebee field must be constant through V'(B 2 ± 6 2 ) = 
(3£/2ft)i7g (if B = 0 this relation does not hold, but 
the bumblebee remains constant nevertheless). Using 
B = Bq = const on Eqs. ( 8 ) and (10) yields 

HS = W%) ' <>=-^ = 0 . ( 19 ) 

As in the previous section, this corresponds to a Universe 
without matter and dominated by the bumblebee field 
alone, which acts as a form of dark energy. 

Given a particular potential V, Eq. (18) and the above 
fully determine both Ho and the nonvanishing value of 
B q. Anticipating the next study, one considers a power- 
law potential of the form 

T)V 

V(B 2 ±b 2 ) = M 4 ~ 2n (B 2 ±b 2 ) n -» V 7 = —r—, (20) 

B z dz b z 


where M has dimensions of mass. Using Eqs. (19) and 
(18) for B 7 ^ 0, one obtains 


B 


2 

0 


Vo 




2 n =F £6 2 

f(l + 2 n) 

/ 2 n 1 ± £ 6 2 

\1 + 2n £il7 2 


M 4 -> 


2 UK ( 2 n 1 ± gfr 2 \ n 2 

3 f V 1 • 2n £M 2 ) 


( 21 ) 


If the bumblebee field vanishes, then one trivially obtains 

Hq = nVo/3. 

Notice that in the above one does not simply find that 
B 2 = =p b 2 and V(B 2 ± b 2 ) = U(0) = 0, as could be 
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naively expected if one assumed that the cosmological 
dynamics should evolve the bumblebee field until it rests 
at its nonvanishing VEV < B 2 >= b 2 : indeed, following 
the discussion of the previous section, the dynamical ef¬ 
fect of the coupling term ^B^B 1 'R^ v found in the action 
functional (1) on the field equations can be interpreted as 
a friction term that arrests the evolution of the bumble¬ 
bee field, with the dissipated energy acting to counteract 
the gravitational attraction of matter and drive an accel¬ 
erated expansion of the Universe. 

One can now look back into the on-shell form of the 
action (1), which — with the prescription (4), B(t) = Bq 
and the de Sitter solution (19) — is reduced to 


S = 



R — 6 




y/^gd A x 


,( 22 ) 


We notice that the dimensionless variable X 2 is directly 
proportional to B 2 , so that no expansion is performed 
around the nonvanishing VEV b 2 ^ 0 arising after the 
LSB; conversely, this definition appears to indicate that 
the vector field is expanded around its symmetric phase 
B 2 = 0. However, since the set of fixed points does 
not change with a shift in the definition of £ 2 , there is 
no distinguishable physics in adopting either form (this 
was checked directly by instead using the definition x 2 = 

as 2 - 6 2 )). 

Furthermore, the calculations below reveal that no 
fixed point arises with B 2 = b 2 (i.e. X 2 = fb 2 ): the 
bumblebee field never rolls to its nonvanishing VEV, due 
to the effect of the coupling with the Hubble parameter 
(and its derivative) in the equations of motion. 

The modified Friedmann equation (8) reduces to 


which, as noted in the preceding section, corresponds 
to the Einstein-Hilbert action with an added constant 
term which acts as a cosmological constant A: if the 
bumblebee field vanishes, the latter is given by its usual 
definition A = 3i£g, by comparison with the usual La- 
grangian density C = (R — 2A)/(2k); if one considers 
/ 0, then it also includes a contribution from the 
coupling between the bumblebee field and the Ricci ten¬ 
sor, A* = 3Hq(1 — ££>o/2) (see Ref. [26] for a study of 
the contribution to a phase of accelerated expansion of 
the Universe arising from a nonminimal coupling between 
matter and the Ricci scalar). 


IV. ANALYSIS OF THE DYNAMICAL SYSTEM 

In order to explore the possible solutions to Eqs. (8-10) 
and confirm the results obtained in the previous section, 
one begins by defining the following dimensionless vari¬ 
ables (see Refs. [27] for a similar treatment in the context 
of quintessence), 


Xi 


kV 

3H 2 


x 2 = £B 2 


= 


fBB 

IT 


(23) 


together with the usual relative matter density and de¬ 
celeration parameter 


n 


M 


Kp 

3 H 2 


q = 


aa 


(24) 


(25) 


Using Eqs. (6), (8) and (10), one obtains 

x[ = 2 (1 + ax 3 + q)x 1 , 
x' 2 = 2x 3 , 

x 3 = (1 - 2q)(l - x 2 ) - 3xi + x 3 (q - 3) , 


where the prime denotes differentiation with respect to 
the number of e-folds N = log a and one defines 


“ fe) S Ffel 




aH 3 


(26) 


1 — x\ + X 2 + xs T , (27) 

while the bumblebee equation of motion becomes 

(2axi + q)x 2 = 0 . (28) 

The two equations above are algebraic constraints that 
must be obeyed by any solution to the dynamical system 
(25). 

The form of Eq. (28) indicates two possible branches, 
corresponding to whether X 2 is vanishing or not. If the 
latter is true, this equation collapses to q — — 2ax\, which 
can then be used to further simplify the dynamical sys¬ 
tem (25). The opposite case X 2 = 0 is more troublesome, 
as it implies that the dynamical system is no longer au¬ 
tonomous [in the sense that x[ = f(xj ) alone]: one would 
have to promote q and to dynamical variables and 
consider the two additional differential equations 

n' M = (2q-l)tt M +/3x 2 -2axi(x 2 +x 3 ) , (29) 
q' = ~/3 + q + 2q 2 . 

Furthermore, it can be shown that the determination of 
the nature of any fixed points would be problematic, as 
the linearization of the Jacobian matrix of the system 
around X 2 = 0 diverges. 

This caveat, however, may be circumvented if one no¬ 
tices that X 2 = 0 is not dynamically relevant if x f 2 ^ 0, as 
the variable just rolls out of its vanishing value and the 
relation q = —2axi becomes immediately valid, allowing 
the aforementioned substitution into Eqs. (25). 

Furthermore, since one is not interested in modeling 
the spontaneous LSB, it is natural to assume that the 
bumblebee field has had time to roll from its previous 
vanishing VEV (before the LSB) towards the broken 
phase < B 2 >= b 2 . 

Thus, one is left only with the pathological case X 2 = 
x' 2 = 0, which is physically hard to interpret, as it would 
imply that the spontaneous LSB has no dynamical effect 
on the Bumbleebee field, which forever rests at its pre- 
LSB vanishing VEV. 

Fortunately, it can be solved analytically; since the 
physical implication of this is that the bumblebee field 
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only imprints an effect on the dynamics through the con¬ 
stant value of the potential Vo, one expects to recover 
the usual picture found in GR for a Universe composed 
of matter and a cosmological constant, where the former 
becomes more and more diluted and the expansion grows 
exponentially at late times. 

Indeed, taking Eqs. (25) and (29) and replacing x 2 = 
x ' 2 = 0 immediately yields X 3 =0, so that 

x[ = 2(l + q)xi , (30) 

3xi = 1 — 2 q , 

= i^q — 1)Om • 

Notice that the last equation is equivalent to the first, 
as can be seen from the Friedmann equation (27), which 
reads = 1 — x% —>• Ll[ m = —x[. One obtains the single 
differential equation 

x[ = 3 (1 — xi) x\ , (31) 


with solution 
xi(N) 
ttm(N) 
q(N) 


1 

1 + Ce- 3JV 

!_ X1= g 

1 C + e 3iv ’ 

l-3 x i_ 3 C 

2 + 2 C + e 3iv 


(32) 


where C is an integration constant. As expected, one 
finds that becomes vanishingly small as the Uni¬ 
verse expands and approaches a de Sitter phase q = — 1 . 
One concludes that, even if the bumblebee field becomes 
locked in the symmetric phase < B 2 >= 0 so that 
X 2 = x ' 2 = 0 , this does not give rise to any unphysi¬ 
cal dynamics, and can proceed towards the study of the 
more relevant scenario X 2 7 ^ 0 . 

As discussed previously, a nonvanishing bumblebee 
field implies that q = — 2 ax\\ replacing this into the dy¬ 
namical system (25) yields the closed set 


The presence of this multiplicative term on two of the 
three dynamical equations (33) renders its behavior very 
important in understanding the evolution of the overall 
system. Since Eq. (27) implies that, for the matter- 
dominated Universe = 1 —>> aq + #2 + #3 = 0, if 
> £b 2 > 0 , x\ will be positive by definition, as seen 
in Eq. (35). In that case, x% will have to be negative, 
causing X 2 to decrease until it crosses the value of £b 2 and 
a diverges (and with it, the trajectories in phase space). 

Therefore, one must enforce X 2 % < £5 2 , where the for¬ 
mer is the initial value of the dimensionless variable x 2 : 
notice that, since one is studying the cosmological dy¬ 
namics after the spontaneous LSB has occurred and the 
potential V has settled into the form (35), this initial 
value is endowed with physical meaning, as it is related 
to the value acquired by the bumblebee field after the 
LSB (which is not modeled directly). 

The inequality X 2 % < £5 2 will imply that a is always 
negative, in order to avoid the divergence of the system: 
as the value of X 2 gets closer to £ 6 2 , a becomes larger 
and (as will be seen numerically), the evolution in phase- 
space becomes more abrupt; thus, it is expected that 
some distance must be allowed between the initial val¬ 
ues of X 2 and £ 6 2 for sufficiently smooth trajectories in 
phase space. Additionally, since the value of £ 6 2 is not 
expected to be very large (since £ is the constant respon¬ 
sible for the Lorentz-violating curvature coupling and no 
evidence of such effects have yet been observed), the im¬ 
posed restraint on the initial values X 2 < £5 2 becomes 
very demanding. 

Under this potential, the system exhibits the fixed 
points shown on I: in particular, fixed points x# and 
xc correspond to the de Sitter solutions found in section 
III B — the former corresponds to either a vanishing cou¬ 
pling £ or a bumblebee field frozen at its unbroken VEV 
B 2 = 0, and is thus an unnatural candidate for a de Sitter 
attractor; conversely, xc displays an additional contribu¬ 
tion due to the coupling between the bumblebee field and 
the Ricci tensor. 


x[ = 2 (1 + axs — 2axi) x 1 , (33) 

x 2 = 2 x 3 , 

x ' 3 = (1 + 4axi)(l — X 2 ) — 3x\ — x 3 (S + 2 axi) , 
and the relation 


(Xl ,X 2 ,X 3 ,Um,<?) _ 

XA (0,0,0, 1, |) 

x B (1,0,0,0,-1) 

Xr , ( 2n-£b 2 q q _A 

y l+2n ’ 1 + 2 n ’ U ’ U ’ L ) 

X D (0,1, 0,0,0) 


X\ + X 2 + 3?3 + UjVf — 1 • 


(34) TABLE I: Fixed points of the model, Eq. (33). 


Following the results obtained in the previous section, 
one adopts a power-law potential of the form ( 20 ), here 
rewritten as 


V{x 2 ) = M${x 2 -(ib 2 ) n , (35) 


where the sign is fixed and £ 6 2 is allowed to have negative 
values. It is worthy of note that, for this potential, the 
variable a becomes: 


n 

x 2 - £6 2 


(36) 


The fixed point x^ corresponds to a matter dominated 
Universe before the bumblebee-induced dark energy be¬ 
comes relevant. This transition requires that this fixed 
point is repulsive, with X# and/or xc attractive. 

The fourth fixed point, X£>, is unphysical, as it cor¬ 
responds to a Universe without matter and only a con¬ 
stant bumblebee field, which either expands or contracts 
linearly or is static (so that q = 0 ); although it is not 
mandatory for this fixed point to be repulsive (as it suf¬ 
fices that our Universe undergoes a trajectory in phase 
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space sufficiently far from it), it is a desirable feature. 

The analysis of the stability of the system can be 
achieved from studying the eigenvalues of the Jacobian 
matrices pertaining to the linearised system in each fixed 
point, which are: 


9 _q _ 4n 4 n 

U ° £b 2 £6 2 

0 0-1 0 

0 2-3 0 

0 0 0 -1 


(37) 


Denoting the three roots of each polynomial tb, i, re, 2 , 
vb ,3 and rc, 1 , re, 2 , re ,3 respectively, the eigenvalues for 
Db and Dc become: 


A b = 


4 n r B , 1 ^,2 A 

& ’ m 2 ’m 2 ’m 2 ) 


(45) 


and 


A c 



rc,i __ r c , 3 \ 

n (£ 6 2 — 1 ) ’ n (£ b 2 — 1 ) ’ n (£ 6 2 — 1 )) 


(46) 


The constraints on n and £b 2 imposed by the stability of 



The eigenvalues \ for matrices Da and Db do not de¬ 
pend on n or £ 6 2 : 

A a = (- 2 ,- 1 ,- 1 , 2 ) , (41) 

A d = (-2,2,-1,-1) , (42) 


- 1.0 - 0.5 0.0 0.5 1.0 

FIG. 1: Parameter constraints assuming x# and xp are 
unstable (white), only xp is stable (light grey) and both 
X 5 and xc are stable (dark grey). The points used in 
the analysis of phase-space trajectories are highlighted. 


and, since at least one of their components has a positive 
real part, the two associated fixed points are unstable. 
For the remaining two fixed points, the issue of stability 
will vary, depending on the values taken by n or £ 6 2 . The 
specific form of the dependence is quite complex, but it 
can be understood in a fairly simple way by defining the 
third degree polynomials (of variable A): 


PbW = A 3 + £ 6 2 (£& 2 - lOn) A 2 + (43) 

2 p6 2 ) 2 [4 n + 4 n 2 - 17<6 2 - 2 (& 2 ) 2 ] A - 


4(£& 2 


8n 


2 n£ 6 2 — 16n 2 £6 2 + (£& 2 ) J 


pc( A) = A 3 + 12n 3 (£& 2 - l ) 3 + (44) 

2 n (£b 2 - 1) (3A + £ 6 2 - l) A + 

3n 2 (^ 6 2 — l ) 2 (5A + 2^6 2 — 2) . 


It can be seen immediately that only the region — 1/2 < 
n < 0 warrants the instability of all of the system’s fixed 
points, regardless of the value of ^b 2 . Additionally, if n 
and £Jd 2 have the same sign, then only one stable fixed 
point is to be found (xp); this fact being especially rel¬ 
evant since the simpler and more readily interpretable 
cases are those for which n is positive (and small). 

If one allows £ b 2 to assume negative values while n re¬ 
mains positive, then the two possible cases of only xp 
being stable or both x# and xp; being stable are sep¬ 
arated by a line that can be roughly described by the 
equation £6 2 ~ —4.76n, with the former lying below the 
line and latter, above. What this means is that to ensure 
the stability of x#, n must increase with £ 6 2 . A similar 
behavior can be noted if, on the other hand, n assumes 
negative values (lower than — 1 / 2 ) and £ 6 2 is positive, 
although in this case the equation describing the divid¬ 
ing line appears shifted (closer to £6 2 ~ —4.76n — 2.38) 
and the area in which both fixed points are stable only 
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asymptotically approaches the axis at £6 2 = 0. Again, 
this is worthy of note since, as mentioned before, £6 2 is 
expected to be small, £ b 2 1. 


V. TRAJECTORY AND EVOLUTION 
ANALYSIS 

Following the previous results, one now studies the tra¬ 
jectories of the variables in the phase-space, especially 
VLm and q , i.e. the bumblebee field, matter density and 
deceleration parameter, respectively - which better pro¬ 
vide a physical picture of the conditions of the Universe 
at the given time. One thus numerically integrates the 
system (33), with initial values corresponding to the de¬ 
parture from the matter dominated era, ~ 1 and 
q ~ 0.5 (the evolution is not affected by choosing ini¬ 
tial conditions sufficiently close to these values). A note 
should be made concerning the initial value of or X 2 i- 
The default would be to start the numerical simulation 
at X 2 = 0, but since the bumblebee field is in that posi¬ 
tion only before symmetry is broken, it is reasonable to 
expect a deviation from the origin, albeit small. 

A. Example of system convergence 

The case where n = 4 is depicted in Figs. 2-4, showing 
that the system collapses entirely to the single attractor, 
as determined by the constraints displayed in Fig. (1). 



FIG. 2: Projection of the phase-space into the &M-q 
plane for n = 4 (two darker plots) and n = 6 (two 
lighter plots). Fixed points are marked as dots in dark 
grey (n = 4) and light grey (n = 6). For both cases 
£6 2 = 10 -2 and X 2 % = 10 -3 (darker), 10 -12 (lighter). 

These trajectories describe a Universe that transitions 
from an initial stage of matter domination (Q,m ~ 1, q ~ 
0.5) to a dark energy dominated era (Om ~ 0, q ~ — 1). It 
is interesting to note that the behavior indicates that — 
instead of evolving smoothly and monotonically towards 
a de Sitter phase — the system overshoots that mark and 
crosses below q < — 1, before asymptotically approaching 
q = — 1. Since the current value of Om ~ 0.3 is well 



FIG. 3: Projection of the phase-space into the X 2 -q 
plane for n = 4 (two darker plots) and n = 6 (two 
lighter plots). Fixed points are marked as dots in dark 
grey (n = 4) and light grey (n = 6). For both cases 
£6 2 = 10 -2 and X 21 = 10 -3 (darker), 10 -12 (lighter). 

before that overshoot and cosmographic surveys indicate 
that q has been steadily declining [28], it is reasonable to 
assume that the lowest value of q lies still in the future. 

Furthermore, because of the discussion surrounding 
Eq. (36), it can be seen that the restrictions imposed 
on the initial values of X 2 drive this overshoot deeply 
into negative values, q <C — 1: mathematically, this is¬ 
sue could be settled by allowing X 2 < 0, which allows for 
much smoother trajectories (if still overshooting q = — 1), 
but results on an added strain on the restriction Eq. (27), 
causing VLm to dip into unphysical (negative) values. 



FIG. 4: Projection of the phase-space into the X 2 
plane for n = 4 (two darker plots) and n = 6 (two 
lighter plots). Fixed points are marked as dots in dark 
grey (n = 4) and light grey (n = 6). For both cases 
£6 2 = 10 -2 and x 2 i = 10 -3 (darker), 10 -12 (lighter). 

Fig. 4 shows that the trajectories decrease mono¬ 
tonically with a physically meaningful, positive matter- 
density profile. Two linear phases are discernible: a very 
steep phase corresponding to the initial plunge of the de¬ 
celeration parameter and consequent overshoot q < 1, 
and a less pronounced decrease during the compensation 
of that overshoot and evolution towards the final state 












X 2 — (1 + £5 2 )/(1 + 2 n). 

B. Example of system divergence 

The second case, for n = 2, is depicted in Figs. 5, 6 
and 7, and exhibits a divergence of the trajectories. In 
this case, the system’s attractor is insufficient to coun¬ 
teract the repulsive effect of the repulsor close to the 
system’s initial values at the start of the simulation. 



°r 


-1 Lf_X_X_X_X_I_X_X_X_X_I_X_J 

0.0 0.5 1.0 

Q-m 

FIG. 5: Projection of the phase-space into the ft m~Q 
plane for n = 4 (two darker plots) and n = 6 (two 
lighter plots). Fixed points are marked as dots in dark 
grey (n = 4) and light grey (n = 6). For both cases 
£b 2 = 10 -2 and X 2 i = 10 -3 (darker), 10 -12 (lighter). 
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q l 

o : 

-l Lf_ji_i_j_i_i_ j _ i _s_j_ l 

0.0 0.5 1.0 

FIG. 6: Projection of the phase-space into the X 2 -q 
plane for n = 4 (two darker plots) and n = 6 (two 
lighter plots). Fixed points are marked as dots in dark 
grey (n = 4) and light grey (n = 6). For both cases 
£b 2 = 10 -2 and X 2 i = 10 -3 (darker), 10 -12 (lighter). 

In this case the analysis is much more straightforward: 
the fixed point is so repulsive that it causes the tra¬ 
jectories to steer clear of the close vicinity of the basin of 
convergence of the attractor, and therefore no physically 
significant convergence towards the latter is attained, as 


the trajectories lead the variables into physically mean¬ 
ingless values. Even if convergent trajectories for this 
case indeed exist, finding them constitutes largely a prob¬ 
lem of fine-tuning. 



FIG. 7: Projection of the phase-space into the X 2 -&m 
plane for n = 4 (two darker plots) and n = 6 (two 
lighter plots). Fixed points are marked as dots in dark 
grey (n = 4) and light grey (n = 6). For both cases 
^b 2 = 10 -2 and X 2 i = 10 -3 (darker), 10“ 12 (lighter). 


C. Overshoot variation 

As mentioned before, the analysed stable trajectories 
present a behavior for the deceleration parameter q over¬ 
shooting its final convergence value at the respective fixed 
point for the chosen potential. The magnitude of this 
effect, though it may appear to depend on the initial 
value taken for x 2 , indeed flattens out after an initial 
phase (for approximately —5.5 < X 2 i < —2) and re¬ 
sponding linearly, eventually stabilising at a final value 
of QMin ~ —15.78, for n = 4, as can be seen in Fig. 8. 



FIG. 8: Variation of the minimum of q with the initial 
value of X 2 for n = 4. 

This result seems to suggest two zones of effective 
system evolution: the first one [—5.5 < log 10 (^2i) < 
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log 10 (£5 2 )] more dynamical and vulnerable to the initial 
values, and the second one [log 10 (^ 2 i) < —5.5] more sta¬ 
ble and robust to variations of the initial conditions. 

As for the other studied case, n = 6, presented in Fig. 
9, it depicts a similar behavior, with the differences be¬ 
ing a faster stabilization [a linear response is seen for 
log 10 (x 2 i) < —4.2, leveling out at qMin ~ —20.47], com¬ 
pared to when n = 4. 


Min (q) 



FIG. 9: Variation of the minimum of q with the initial 
value of X 2 for n = 6. 


VI. CONCLUSIONS AND OUTLOOK 

In this work, one has considered the bumblebee model, 
a vectorial extension of general relativity under a spon¬ 
taneous symmetry breaking mechanism. Since the study 
is made considering that the bumblebee vector is only 
nontrivial in the time component, which depends only 
on time, no additional breaking of Lorentz symmetry oc¬ 
curs, aside from the one derived from adopting the FRW 
metric, which is equivalent to a foliation of space-time 
into spatial pictures labelled by a cosmological time. 

It was found that the model predicts four points of 
equilibrium, in agreement with the preliminary study of 
the equations; two are unstable regardless of the value of 
the model’s parameters, corresponding to the static and 
matter-dominated cases. The remaining two yield a Uni¬ 
verse undergoing an accelerated expansion, which can be 
ascribed to the arise of a quantity akin to a cosmological 
constant: this can be caused by the effects of the po¬ 
tential alone (for practically all values of the parameters, 
with the exception of —4 < n < 0), making the identifi¬ 
cation A = 3Hq 2 , or to the effects of a concerted evolution 
between the potential and the bumblebee vector’s non¬ 
zero component, resulting in A* = —3Hq 2 (1 -f Vc) /2, 
where Vc = (1 ± £6 2 )/( 1 + 2 n) is the value of the poten¬ 
tial on the fixed point x^. 

It was noted that, when considering a power-law type 
of potential, even within the parametrical constraints 
that guarantee the existence of one attractor, a conver¬ 
gent scenario isn’t always found: in particular, this was 
not attained in the case of odd powers, which, for all 
tested values, appear to cause the system’s variables to 


be strongly repelled initially, thus diverging from the de 
Sitter attractor. This can be explained by the fact that, 
for odd n, the variable x\ will be initially negative if 
X2 < £5 2 and positive if X2 > £,b 2 . But the discussion 
concerning Eq. 36 guarantees the system will diverge 
in the latter case, and in the former, the negativity of 
x\ disturbs the balance of the Friedmann Equation (27), 
disallowing convergence. 

The above does not guarantee that the system does 
not allow converging trajectories, only that they are not 
dense in the (continuous) set of all possible trajectories 
with physically significant starting points, and even if 
one was to be found, there would be no basis to argue 
the great fine-tuning needed for its manifestation. 

Additionally, the analysis of convergent trajectories in 
phase-space predict that q always overshoots the final 
value of —1, but since the present-day Universe is still 
on the descending slope of the trajectory, prior to that 
overshoot, such behavior constitutes a prediction of the 
model, rather than an experimental constraint (albeit one 
not to be verified in the foreseeable future). 

Expectations concerning £6 2 , X 2 and convergence of the 
system may, however, be used in restricting the window 
of accessible values to the dynamical system’s variables: 
the major factors in that process can be categorized into 
three primary groups: (1) regularity of the solutions, (2) 
elimination of unphysical behavior and (3) mathematical 
stability of dynamical system. 

1. n < — 1/2 or n > 0: This bound is required from a 
stability analysis of the fixed points of the dynam¬ 
ical system, as there are no converging trajectories 
for values of the exponent n outside the indicated 
window. However, it can be argued that there is no 
guarantee that the evolution of the Universe must 
lead to stable convergence — or that, even in the 
presence of suitable attractors, a region of converg¬ 
ing trajectories can be found for a specific choice 
of model parameters n and £6 2 without fine-tuning 
the initial conditions. 

2. X 2 i < £5 2 : After the spontaneous LSB has occurred 
and the potential has settled into its form (35), 
this upper bound for the value of the bumblebee 
field must be obeyed, otherwise the system diverges 
through the action of a. 

3. X 2 i < £5 2 : The stronger restriction, not contem¬ 
plated in the Friedmann Equation, is related with 
the behavior of the deceleration parameter g, and 
ascertains that the minimum of the solution for this 
variable (the so-called overshoot) shows a tendency 
to explode if the absolute initial values of X 2 and 
£6 2 are too close. Once more, this is due to the di¬ 
mensionless quantity a diverging and dominating 
the dynamic of the system, making impossible the 
smooth transition required in order for there not to 
be an overshoot. 
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Regarding the obtained interval for the exponent n 
where at least one of the de Sitter fixed points is an at¬ 
tractor (first point of the list above), the comparison with 
a more complex, nonpolynomial potential could provide 
interesting hints: indeed, if the latter can be expressed 
analytically as a series of powers of R 2 , the two should 
match in the relevant regime for B 2 where the latter can 
be truncated. 

Such a potential naturally arises from the asymptoti¬ 
cally free solutions of a quantized bumblebee model [15], 
being expressed in terms of the Kummer (confluent hy¬ 
pergeometric) function M(k — 2,2 ,z) — 1. However, the 
stability analysis found there did not assess the dynam¬ 
ical (cosmological) evolution of the bumblebee field, but 
concentrated on the existence of minima of the potential, 
so that the former could rest on a nonvanishing VEV — 
with excitations around the latter giving rise to Nambu- 
Goldstone modes associated with the ensuing LSB. 

Our scenario differs significantly: it is classical in na¬ 
ture, has one nonvanishing minima defined by the power- 
law form of the potential, and the fixed points show that 
the bumblebee field does not evolve towards this nonva¬ 
nishing VEV; furthermore, the results obtained in this 
study rely extensively on the presence of a coupling with 
geometry, which is absent from Ref. [15]. As such, no 
direct comparison is feasible. 

As a final remark, we note that the fairly rich struc¬ 
ture of the model prompts more focused studies: going 
beyond further inquiries into cosmological dynamics [e.g. 
by selecting a different potential V(B 2 ±b 2 ) or including 
a coupling between the bumblebee field and the Ricci 
scalar], one anticipates that interesting results should 
arise from a spatially oriented bumblebee that breaks the 
homogeneity and isotropy postulated by the cosmologi¬ 
cal principle: this could bear a possible relation with the 
putative “axis of evil” reported in the cosmological back¬ 
ground radiation [29]. 

A different approach could also consider several local 
instances of spontaneous Lorentz symmetry breaking, in 
analogy with Condensed Matter Physics [30]: one can 
argue that this is compatible with the above-mentioned 
possibility of a spontaneous LSB in a spatial direction, so 
that the bumblebee field acquires a particular (random) 
orientation only within a region with a well defined co¬ 
herence length, beyond which another orientation is ran¬ 
domly adopted during the LSB. If the size of the observ¬ 
able Universe is much larger than the coherence length, 
it is physically plausible that the large scale superposi¬ 
tion of several regions with different spatial orientations 
for the bumblebee field gives rise to a globally homoge¬ 
neous and isotropic Universe. The collision of different 
“bubbles” could also produce interesting physics and hy¬ 
pothetically lead to observable relics of the LSB. 

The in-depth study of such a mechanism would most 
probably require modeling the actual LSB (i.e. the evolu¬ 
tion of the shape of the potential, so that the VEV of the 
bumblebee field becomes nonvanishing) and encompass 
similar considerations to the well-known Kibble-Zurek 


mechanism [30], albeit with the added complexity of con¬ 
sidering a vector field instead of a scalar. By the same 
token, one remarks that the interface between regions 
with different spatial LSB could give rise to topologi¬ 
cal defects, which would physically correspond to areas 
where Lorentz symmetry holds. 

However mesmerizing, the above scenario clearly re¬ 
quires heavy-duty numerical computations, and is thus 
incompatible with the stated purpose of this work: to 
clearly identify the possibility of driving an accelerated 
expansion of the Universe via a nonminimally coupled 
bumblebee vector field. 
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Additional bumblebee-curvature coupling 


As discussed in section III, the action functional (1) 
considered in this study follows Ref. [9], with only a cou¬ 
pling between the bumblebee field and the Ricci tensor 
considered: we argue that this captures the essential fea¬ 
tures of the model, namely the ensuing de Sitter attractor 
xc depending on the coupling strength £. 

If one additionally considers the coupling term in the 
action functional 

s = j ^B^R^d 4 x , (47) 

where x is a coupling strength (with dimensions [x] = 
[£] = M -2 ), the modified bumblebee Eq. (3) becomes 

= 2 (v'B u - _ 2L b v r}j , (48) 

and the r.h.s. of the modified Einstein field equations 
acquires the extra terms 


—x[BaB a G /Jb y + RB^B V + (49) 

g^n(B a B a )-x/^y(B a B a )} . 

Introducing the FRW metric (5), Eq. (48) becomes 



3(2x + Q a _ 3X H 2 

2k a k 


B = 0 , 


(50) 


so that, considering the additional contribution to the 
derivative of the potential the curvature appearing on 
the r.h.s. of the field equations, AV' = (x/2ft)R, the 
scalar curvature transforms as 


R^(l-ZB 2 )R-3 X n(B a B a ) , (51) 
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The modified Friedmann and Raychaudhuri Eqs. (8,9) 
thus read 

H 2 [ 1 - (£ + x)B 2 ] = \k{ P + V) + (C + 2 X )HBB , 

( h 2 + 2 1 ) [i-(x+e^ 2 ] = 

/W + 4(£ + x)-£™ + (£ + 2 X ) (B 2 + BB) . (52) 


tween the bumblebee field and geometry. The modified 
bumblebee Eq. (50) becomes 

v*- «* + ° * » ( 53 ) 

again providing an intrinsic relation setting B{t) = Bq = 
const.. Introducing this relation into Eq. (52) yields 


One cannot absorb the effect of the additional coupling 
into a redefined parameter, as both factors £ + x and 
£ + 2x appear; notwithstanding, the above shows that the 
dynamics do not differ significantly from those obtained 
with x = 0 — as no additional combinations of H, B or 
its derivatives arise. 

One can apply the procedure followed in Section III B 
to ascertain the possibility of obtaining a de Sitter phase 
from this more evolved scenario where two couplings be¬ 




kV 0 


3[1 -« + x)Bg] 


P = 0 , 


(54) 


which, as before, describes a Universe with no matter 
content and with the potential term acting as a dark 
energy component driving its accelerated expansion. In¬ 
terestingly, the above shows that in the particular case 
£ = —x the two couplings cancel. 
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